Antiferromagnetism and charged vortices in high- Tc superconductors 
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The effect of the long-range Coulomb interaction on charge accumulation in antiferromagnetic 
vortices in high- Tc superconductors is studied within a Bogoliubov-de Gennes mean-field model of 
competing antiferromagnetic and d-wave superconducting orders. Antiferromagnetism is found to be 
associated with an accumulation of charge in the vortex core, even in the presence of the long-range 
Coulomb interaction. The manifestation of ll-triplet pairing in the presence of coexisting d-wave 
superconductivity and antiferromagnetic order, and the intriguing appearance of one-dimensional 
stripelike ordering are discussed. The local density of states in the vortex core is calculated and is 
found to be in excellent qualitative agreement with experimental data. 
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I. INTRODUCTION 

The vortex cores of high-T^ superconductors have been 
the subject of intense investigation in recent years. Ad- 
vances in experimental techniques, especially in the area 
of scanning tunneling microscopy (STM)fi4242i4 have al- 
lowed unprecedented exploration of the electronic states 
at the atomic scale. Local information on vortex core 
states has also been provided by high-field nuclear mag- 
netic resonance (NMR) experiments^S'SiS'S which use in- 
homogeneity in the magnetic field to extract a spatially 
resolved relaxation rate. 

Among the principal results of these experiments is the 
suppression of the low-energy density of states in the vor- 
tex corejSi with either an absence, or possibly a strong 
splitting of the zero-bias conductance peak (ZBCP)fi'2ii 
indicating some form of ordering in the vortex core. STM 
results also show the presence of low-energy vortex core 
states which extend beyond the vortex core radiusfi*i4 
as defined by the superconducting coherence length ^. 
These extended vortex-induced states were further ob- 
served to form a localized checkerboard pattern along 
the CuO bonds with a period near 4.3ao (Ref.i). A sim- 
ilar checkerboard pattern with a period close to 4.5 oq 
has also been seen in more recent STM experiments in 
zero field, both in the pseudogap stateiS and in the su- 
perconducting stateiiii^of BiaSrsCaCuaOg+a (BSCCO). 
These observations are corroborated by STM experi- 
ments on the superconductor NaajCaa-xCuOaCla (Na- 
CCOC) which show strong ordering with a period of 4ao 
in a pseudogaplike state. It has been proposed that this 
checkerboard pattern is due to the formation of a pair- 
density wave of phase-incoherent Cooper pairsj ^'^d^d^ 
however, this interpretation remains controversial and 
the existence of other orderings, such as dynamically 
fiuctuating stripes, spin-density wave orderingirS or 
a particle-hole charge-density wave^^S remains a likeli- 
hood. Indeed, neutron-scattering experiments have re- 
ported the existence of a periodic modulation of antifer- 
romagnetic correlations in zero fieldj^-^'^^i^^i^'^i^^i^^i^'''i^^i^^ 

Many experiments have also reported observations of 



field-induced magnetic order in the high-Tf, supercon- 
ductors. It has been noted that the periodicity of the 
charge ordering observed in STM (Ref. 0) could be con- 
sistent with inelastic neutron-scattering experiments on 
optimally doped La2-£!;Sra;Cu04 {x = 0.163) that show a 
field-induced signal in the low-frequency spin-fiuctuation 
spectrum, near the (tt, tt) point in reciprocal space, sug- 
gesting the existence of a fluctuating spin-density wave 
with a periodicity close to Suq along the CuO bond direc- 
tions and a correlation length i > 20 lattice spacingsi^ 
More local measurements using spatially resolved NMR 
also provide strong evidence for the presence of antifer- 
romagnetic order in the vortex cores of near-optimally 
doped YBazCuaOy.a (Ref. H) and TlsBasCuOg+a (Ref. 
0). Muon-spin resonance measurements of the mag- 
netic field distribution of underdoped YBaaCuaOg.so re- 
veal structure in the high-field tail that is not seen in 
optimally doped YBaaCusOg.gs and is consistent with 
static antiferromagnetism."^^ Elastic neutron-scattering 
has shown a strong enhancement of incommensurate 
static antiferromagnetism in an applied magnetic field in 
underdoped La2-a;Sr^Cu04 with x — 0.10 (Ref. Is^ and 
with X = 0.12 (Ref.l^. The antiferromagnetism (AFM) 
is observed to have a very long correlation length > 400 
A showing that AFM and superconductivity (SC) coexist 
in the bulk."^^ Similar results are seen in elastic neutron- 
scattering experiments on La2Cu04_|_y (Ref. Isl) . The 
picture that emerges from these various measurements 
is that static or dynamic antiferromagnetic order is in- 
duced, or at least strongly enhanced, by an applied mag- 
netic field. This antiferromagnetic order appears to be 
nucleated at the vortex cores, and extends into the bulk 
of the superconductor with a correlation length that is 
much longer than the superconducting coherence length. 

Charge-density modulations induced at vortex cores 
by the presence of antiferromagnetic or spin-density- 
wave order have been predicted within various theoretical 
approaches ii2iSi^Si2L^2i22ii2iii In partial agreement with 
the experimental results discussed above, common fea- 
tures of many of these calculations include an antiferro- 
magnetic or spin-density-wave order that is nucleated at 
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the vortex core (where d-wave superconductivity is most 
strongly suppressed) and which extends well beyond the 
core region into the bulk of the superconductor. This an- 
tiferromagnetic order is accompanied by a local modula- 
tion in the charge density that reaches a maximum at the 
centre of the vortex core and is seen to approach half fill- 
ing in some cases. ■^^■^'^ The vortex core charge is predicted 
to have a strong doping dependence!^ In the overdoped 
region of the phase diagram the vortex cores are not an- 
tiferromagnetic, and the difference between the chemical 
potentials of the normal core and the surrounding super- 
conducting region leads to a small decrease in the charge 
density at the corei^Si^ In this case the vortices are pos- 
itively charged. As the doping is decreased, the vortex 
cores become antiferromagnetic and expel holes into the 
bulk of the superconductor, causing the vortices to be- 
come negatively charged. 

Since there is no reason to expect pair-density-wave 
ordering to push the local charge density toward half fill- 
ing, the presence of negatively charged vortices in the 
underdoped cuprates could be a distinguishing feature of 
field-induced antiferromagnetic ordering. However, the 
ability of this induced charge to survive the effects of 
the long-range Coulomb (LRC) interaction has not been 
established. In order to determine whether negatively 
charged vortices can persist as a robust prediction of 
a theory containing antiferromagnetically ordered vor- 
tices, we have set out to include the long-range Coulomb 
interaction in a self-consistent mean-field theory of the 
cuprates that contains both antiferromagnetic and su- 
perconducting order parameters. The self-consistent ef- 
fect of the long-range Coulomb interaction on the doping 
and magnetic field dependence of the antiferromagnetic 
order and vortex-core charge is calculated. Our results 
suggest that an antiferromagnetic vortex core is always 
associated with an accumulation of electrons in the core 
region, even in the presence of the long-range Coulomb 
interaction. However, both the AFM order and the as- 
sociated vortex core charge are weakened with increased 
LRC interaction strength, and eventually vanish simul- 
taneously above a critical value. We further find that the 
local density of states (LDOS) in the vortex core is modi- 
fied by the long-range Coulomb interaction in such a way 
that the essential features of STM data on YBCO (Ref. 
^ are captured within a simplified model of competing 
AFM and dSC orders. We also study the manifestation 
of H-triplet pairing in the presence of coexisting dSG and 
AFM order, and discuss the intriguing appearance of one- 
dimensional stripelike ordering. 



II. MODEL AND METHOD 

We describe a two-dimensional d-wave superconductor 
in an external magnetic field by the t-J Hamiltonian with 



an additional long-range Coulomb interaction 

where V — e^/e/jO^ with dielectric constant eh, and \rij\ 
is in units of the lattice constant oq. The LRC inter- 
action is added to the Hamiltonian in order to study 
the competing effects of the local antiferromagnetic or- 
der in the vortex cores and the Coulomb repulsion be- 
tween electrons J -i ^ -' i' We start with the mean-field model 
of Ghosal et al.^^ and include the effects of the long- 
range Coulomb interaction over all sites for the Hartree 
shift (which includes the density-density interaction) and 
at the nearest-neighbor level for all other order param- 
eters. It is important to note that screening effects 
will be accounted for by self-consistent calculation of 
the Hartree shift — which allows for rearrangement of the 
charge density — as in density functional theory (a simi- 
lar approach has been taken in a zero-field study of the 
t-J model by Arrigoni et alM). With the LRC interac- 
tion incorporated in this way, the effective Hamiltonian 
of Ref. 39 becomes 




+ E {^A^4+s^-. + H.c.) , (2) 

where rf ^ is the Fock shift that renormalizes the hopping 
amplitude t (we set t = 1 in this paper), and rji is the 
Hartree shift that renormalizes the chemical potential /z. 
Here S = ±x, ±y, the four nearest neighbors of any site, 
and (l)f = (tt/^o) J^^ Mr) ■ dr. The local dSC (Af ,,) 
and AFM {m.i) order parameters satisfy the following 
self-consistency relations: 




(3) 

and 

= (4) 

where Q = (tt, tt) is the antiferromagnetic wave vector. 
The Fock shift 

(5) 

is assumed to have no dependence on the direction of hop- 
ping across a given bond, and is allowed to renormalize 
only the magnitude and not the phase of the hopping en- 
ergy t. This is consistent with the assumption employed 
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here that the magnetic field is uniform throughout the 
sample. The Hartree shift 

where the sum in j is over all sites in a magnetic unit 
cell, contains the density-density part of the long-range 
Coulomb interaction. The constants aij are calculated 
by the method of Ewald summation with the assumption 
that the magnetic unit cell is periodically repeated in all 
three directions in space, separated by a distance in the 
perpendicular direction of c ~ 3ao, and surrounded by 
a uniform, neutralizing background charge with a static 
dielectric constant e^. In this paper, we report on re- 
sults for V from to 0.35. For t ^ 0.5 eV, this cor- 
responds to a static dielectric constant for the bound 
charges from en ^ 00 down to eh ~ 20. Since we 
are working at high fields, we assume a square vortex 
lattice^^*^ with nearest-neighbor vortices along the Cu- 
O bond directions: our results, which focus on the vortex 
core structure, do not depend on this choice. 

In order to diagonalize the effective Hamiltonian, we 
make a Bogoliubov-de Gennes (BdG) transformation 
and block diagonalize the resulting BdG Hamiltonian 
by exploiting the magnetic-translational symmetry of the 
problem. Starting with an initial guess for all local vari- 
ables and for the value of /i, we solve for all eigenvalues 
and eigenvectors of the BdG matrix for each magnetic 
wave vector k. Results in this paper are reported for 
24 X 12 wave vectors. Note that the time-reversal sym- 
metry of the BdG Hamiltonian,^'' E„{\i) — k), 
allows us to find the eigenvalues and eigenvectors of the 
k-dependent blocks of the BdG matrix in (k, — k) pairs 

^,_.(k); 1, <,(-k) ; ■ 

The eigenvalues and eigenvectors are used to recalculate 
all of the local variables using the self-consistency equa- 
tions defined above [Eqs. ©"©I- This iterative process 
is repeated until the largest change in all site-dependent 
variables and in the average charge density (which is con- 
trolled by tuning the chemical potential /i) is less than 
10^^. The addition of the long-range Coulomb interac- 
tion destabilizes this iterative procedure, and we found 
it necessary to develop a modified Broyden's method to 
achieve convergence in the Hartree shift. The details of 
the convergence method are discussed in appendix 1X1 

III. RESULTS 

A. Effect of the long-range Coulomb interaction on 
the vortex core charge 

In the absence of the long-range Coulomb interaction 
the suppression of d-wave superconductivity at the vortex 
cores leads to the emergence of antiferromagnetic order 



with a correlation length that is much longer than 
the superconducting coherence length ^, for the range of 
parameters that is appropriate to the cuprates. The an- 
tiferromagnetism is strongest at the center of the vortex 
core, and creates a non uniform charge density through 
the expulsion of holes from the vortex core into the bulk 
of the superconductor. In this paper, we focus our at- 
tention on the added effects of the long-range Coulomb 
interaction on this charge ordering. 

As one would expect, the principal effect of including 
the Coulomb repulsion is to reduce the magnitude of the 
peak charge density at the center of the vortex core. As 
can be seen in Fig.^ which plots the non uniform charge 
density 5{ni) — (rii) — nave, the large non-uniform elec- 
tron density in the vortex core is first reduced and then 
eliminated as the LRC strength is increased from ^ = 
to V = 0.35. Initially, the effect of the LRC interaction 
is to sharpen the charge density profile and to induce a 
screening region of reduced charge density immediately 
around the vortex core. This screening region separates 
the charge density profile into two distinct regions: inside 
the core AFM order governs the charge density, pushing 
it toward half filling; outside the core AFM no longer 
controls the charge density and there is a weak dip in 
the charge density along the nodal directions beyond the 
screening region, due to the effects of dSC order. 

Interestingly, the shape of the antiferromagnetic order- 
ing (see Fig. [2J) is not affected by these dramatic changes 
in the short-range structure of the charge density. This is 
refiectivc of the long correlation length of the AFM order. 
However, the magnitude of the AFM order is reduced 
along with that of the non uniform charge density. As the 
interaction strength is further increased to ^ = 0.35, the 
charge density modulations are greatly weakened, with a 
small amplitude on the order of 10"'^ electrons per site. 
At this point, the AFM order at the vortex core has been 
completely destroyed (|m| ~ 10"'*) by the reduction in 
the local charge density. In the absence of AFM order, 
the vortex core becomes weakly positively charged due to 
the gap variation between the "normal" core and the sur- 
rounding superconductor The charge density profile 
at V — 0.35 is equivalent to that seen in more highly 
doped systems where AFM order is absent even a.t V ~ 
(discussed below in Sec. IIII CJ|) . 

The d-wave superconductivity (dSC) is also suppressed 
by the Coulomb repulsion, which weakens the nearest- 
neighbor pairing attraction [as can be seen in the second 
part of Eq. ©]. Fig. 121 shows the evolution of the dSC 
order as the Coulomb interaction is increased from V = 
to V ~ 0.35. The overall shape of the dSC order does 
not change; but, the coherence length gradually increases 
as the magnitude is reduced. Note that it is the site- 
dependent gauge-invariant gap operator, defined as 

S,(T 

that is shown. 
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It is clear that the survival of antiferromagnetic order 
at the vortex core is dependent upon the existence of a 
local charge density that is closer to half filling than the 
average density. This dependence is evident in Fig. 01 
which shows the effects of the long-range Coulomb inter- 
action on the strength of the dSC, on the AFM order, 
and on the peak charge density in the vortex core. The 
AFM decreases gradually until V > 0.2, where it be- 
gins to rapidly fall off. Both the antiferromagnetism and 
the charged vortex core disappear between V = 0.3 and 
V — 0.35. The dSC order |Ao| decreases linearly at first 
and then gradually starts to fall off more quickly with 
V, which is consistent with the reduction of the pairing 
interaction strength in the second part of Eq. Q , and is 
able to survive at larger values of V where the AFM and 
the charge order are destroyed. The existence of antifer- 
romagnetic order clearly implies, and is dependent upon, 
the existence of negatively charged vortex cores. 



B. Magnetic field dependence of antiferromagnetic 
order and core charge density 

In an attempt to better characterize the effect of the 
long-range Coulomb interaction on the vortex core charge 
we have also studied the dependence of the AFM order, 
and its accompanying charge density order, on the size 
of the unit cell {H = 2^o/N^Nyal « 34 T for a 20 x 40 
unit cell). At V = 0, the AFM increases and the charge 
density moves closer to half filling with increased unit cell 
size. In a larger unit cell the AFM order, which has a 
correlation length im >> is able to reach larger values 
by spreading out its variations over longer distances. Fig. 
Elshows that the LRC interaction has a dramatic effect on 
this relationship: when V > 0, the AFM and the charge 
density at the center of the vortex core initially increase, 
but then start to decrease with increasing unit-cell size 
above a unit-cell size of 22 x 44. 

As discussed above, the shape of the AFM order does 
not change with increased V, despite dramatic changes 
in the structure of the charge density. The observed re- 
duction in AFM order in larger unit cells is likely to due 
to the fact that the peak in the AFM is spread out over 
large distances. At V = this is to the benefit of the 
AFM; however, when the LRC interaction is turned on 
the charge density peak is sharpened, and a hole rich re- 
gion of screening charge forms around it. As the unit cell 
gets larger, more of this hole-rich region will lie under- 
neath the central peak in the AFM order. We suggest 
that the weakening of the AFM observed in larger unit 
cells is due to this lowering of the 'average' charge density 
under the peak in the AFM order. 

Taken to its limit, this argument suggests that AFM 
will eventually disappear as the magnetic field is low- 
ered, if it maintains the spatial structure commensurate 
with that of the supcrcurrent that is seen in Fig. [2 It 
is, however, more likely that the spatial structure of the 
AFM will change at some point as the unit cell size is 



increased beyond the AFM correlation length pre- 
serving the AFM vortex core. Indeed, such an effect 
is discussed below in Sec. IIII El It is also important 
to note that AFM order with shorter-range 2D spatial 
structure that is not tied to the unit cell size, as seen in 
other ■workfi2i2Si2iiiiii2ii2 could be accessed within our 
model by tuning of the model parameters and introduc- 
ing a next-nearest-neighbor hopping,^^'**^ Such shorter- 
range ordering may well be less perturbed by the screen- 
ing effect of the long-range Coulomb interaction at low 
magnetic fields. 

C. Doping dependence 

We have also studied the dependence of the vortex core 
charge on the average density of electrons, or the dop- 
ing. For V — 0, our results are in agreement with those 
of Chen et al. who observed positively charged vortices 
with negligible AFM order for larger values of the hole 
density, with a transition to negatively charged AFM vor- 
tices as the average density approaches half filling.^" For 
a unit cell of 20 x 40 sites with two vortices this tran- 
sition from positive to negative vortices occurs near an 
average density of Tiavo = 0.855, a,t V — (see Fig. EJ. 
The chief effect of increasing V is to push this transition 
closer to half filling and to reduce the magnitudes of both 
the AFM order and the vortex core charge. The AFM 
region shrinks gradually at first and then rapidly with 
increasing V. This is consistent with the behavior seen 
in Fig. ^ where the AFM at the vortex core is seen to 
decrease more rapidly at larger V . The transition from 
positively to negatively charged vortices occurs when an- 
tiferromagnetism begins to be nucleated in the vortex 
cores. 



D. It-triplet pairing 

It has been shown that Il-triplet orderiSSi^fli^iiS is able 
to develop self-consistently in the presence of coexisting 
dSC and AFM order, even in the absence of any interac- 
tion which generates it directly. The local Il-triplet order 
parameter 

nf = ^{(c.+54Q,t) - (Q,iC.+5,T)}e'«-"-% (9) 

where Q = (tt, tt), has a triplet pairing amplitude in the 
Sz — channel. The prefactor of J/4 is introduced for 
the sake of comparison with the other order parameters. 
We find that the shape of the Il-triplet order does not 
change with increasing V: within the vortex core, jllil 
rises rapidly from zero to its maximum value within a 
superconducting coherence length ^; outside of the vor- 
tex core I Ail is approximately uniform, and the shape of 
jlljl follows that of the AFM order, falling with increas- 
ing distance from the core to zero at the lines of zero 
supercurrent before rising again to smaller peaks in the 
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corners of the unit cell. The magnitude of the Il-triplet 
order decreases with increasing LRC, along with that of 
the AFM order, and falls below lO^^i when V reaches 
0.35. 

The Il-triplet order parameter was introduced by 
ZhangfO in the context of S0(5) theory, which proposes a 
unified theory of AFM and dSG order for the cuprates. It 
was argued on general symmetry grounds that a theory 
containing A.^, 11^, and orders will satisfy the relation 

li*A, = -m,(l - (n,)), (10) 

where the variables without the "tilde" are those of Eqs. 
0, and jnj without the prefactors indicating the 
interaction strength. At = this symmetry relation is 
satisfied within 5% on a site-by-site basis. When ^ > 
we find that Eq. (|10|l is still satisfied within 10-20 % over 
most of the unit cell, except in the vortex core where the 
charge density is pushed away from half filling such that 
(1 — (ricoro)) no longer approaches zero. In the presence 
of the LRC interaction, Il-triplet pairing is observed, but 
the symmetry relation connecting 11^, A^, and the 
local hole density is weakened. 

E. One-dimensional anisotropy 

When the exchange interaction strength J is increased 
above 1.15, the AFM and the charge density sponta- 
neously develop a slight one-dimensional anisotropy at 
V = 0: both the AFM and charge orders extend to 
a greater distance along the x oi y direction (parallel 
to the Cu-0 bonds). We have checked the validity of 
this behavior by repeating the calculations on a square 
40 X 40 site unit cell containing four vortices, and by 
comparing the ground state energies of systems with x- 
and y-oriented anisotropy in both rectangular (two vor- 
tices) and square (four vortices) unit cells. All results 
were found to be equivalent within the limits of numer- 
ical accuracy. Similar behaviour has been seen in other 
BdG studies of the vortex state based on mean-field solu- 
tions of the extended Hubbard model with added nearest- 
neighbor pairing42iS 

We have found that this one-dimensional anisotropy 
is strongly enhanced by the long-range Coulomb interac- 
tion. At V = (not shown) the AFM and the charge 
density orders have an elliptical shape. At ^ = 0.35, 
as shown in Fig. the AFM and the charge density are 
strongly one dimensional, extending throughout the unit 
cell along the y direction while oscillating in the x direc- 
tion with a period of 20ao for the AFM order and lOao for 
the charge density order. The period of this oscillation is 
set by the size of the unit cell. We note that other period- 
icities and more complex orders are accessible within our 
model if terms such as a next-nearest-neighbor hopping 
are addedi^i*^ The dSC order is only slightly affected 
and, for the most part, retains its 2D structure. The 
AFM is strongest in the vortex core and has the opposite 



sign at the edges of the unit cell. The magnitude of the 
AFM at the edges is almost as large as it is in the center 
"stripe" and has, surprisingly, been enhanced by the LRC 
interaction. The charge density tracks the magnitude of 
the AFM and shows none of the rich, short-range, two- 
dimensional structure seen at J — 1.15 in the presence 
of the LRC interaction. With J = 1.3, the main effect 
of the LRC interaction (apart from reducing the mag- 
nitude of the AFM and the charge density at the core) 
is to strongly enhance the one-dimensional ordering that 
was only weakly developed at 1^ = 0. We suggest the 
following reasons for this phenomenology. First, as the 
exchange interaction strength J is increased, the domain 
walls at which the AFM order parameter nii changes 
sign cost increasing amounts of energy. As the AFM 
is strengthened, this energy cost creates a tendency to 
straighten out and shorten these domain walls. We have 
also noticed a much weaker tendency toward anisotropy 
at J = 1.15 in larger unit cells with 24 x 48 and 26 x 52 
sites, that leads to strengthening of the AFM (as can 
be seen in Fig. EJ. This anisotropy is not obvious at 
V = 0, and is only weakly apparent for > 0. This 
unit-cell size dependence suggests that the correlation 
length £rm which decreases with increasing J, may also 
be playing a role. Earlier work on stripes within SO (5) 
theory?- showed that phase-separated AFM and dSC re- 
gions change from a droplet phase to one with alternating 
stripes when the AFM and dSC regions occupy roughly 
equal areas of the superconductor. The anisotropy ob- 
served here may be partly driven by coincidence of £m 
with the length of the unit cell, such that the AFM order 
parameter would like to occupy half of the unit cell area. 
Secondly, the LRC repulsion reduces the charge density 
at the vortex core and creates a hole-rich screening re- 
gion immediately next to the core that screens the vortex 
charge from the rest of the unit cell. The charge that is 
moved away from the core is thus pushed to the edges of 
the unit cell. This balancing of the charge reduces the 
AFM in the center and strengthens it along the edges of 
the unit cell. The LRC interaction does its best to cre- 
ate a uniform charge density and strongly enhances the 
ID character seeded by the AFM order by washing out 
modulations in the charge density along the y direction. 
Further increase of the LRC interaction will maintain this 
stripelike order, but will reduce the magnitude of both 
the charge and AFM oscillations. 



F. Local density of states 

The local density of states provides a basis for com- 
parison of theoretical results with powerful experimental 
tools. As is well known, the zero-bias conductance peak 
(ZBCP) predicted to exist in the vortex core by standard 
BCS theory for d-wave superconductors^ has not been 
experimentally observed. Scanning-tunneling microscopy 
(STM) measurements of the vortex state of YBCO (Ref. 
B) and of BSCCO (Refs. 1311 show a complete absence 
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of the ZBCP in the vortex core. Low-energy structure in 
the conductance at ±5.5 meV in YBCO (Ref. [ij) and at 
±7 meV in BSCCO (Refs. is observed in the vortex 
core spectra, and is seen to persist outside of the core to 
distances much longer than the coherence length. 

We present our results for the local density of states in 
the vortex core region in Fig. |S1 In agreement with pre- 
vious results (^SiSSiSLS we find that the presence of anti- 
ferromagnetic order in the vortex core provides a mech- 
anism for splitting of the ZBCP, both in the presence 
and absence of the LRC interaction. More importantly, 
the LRC interaction changes the spectrum at low energy, 
leading to excellent qualitative agreement with experi- 
mental results. The main features in the local density of 
states attributable to the LRC interaction are the follow- 
ing: the unphysical spin gap at E ^ 0.5 t is washed out 
for V — 0.15 due to weakening of the AFM order; the 
asymmetry in height between the split zero-bias conduc- 
tance peaks is reduced; and the peak splitting is found 
to decrease linearly with |tocoic|, as expected for the ap- 
pearance of an order parameter that breaks time-reversal 
symmetryiSiS2iSiiS2iS^ The bottom trace of Fig. H shows 
the local density of states averaged over all sites within 
the vortex core. Given the simplicity of our minimal 
mean-field model, the agreement at V — 0.15 with STM 
results on YBCO (Ref.[l|) is remarkable. Almost all of the 
qualitative low-energy features of the experimental data, 
including the appearance of slightly asymmetric peaks 
at E ^ A/4, with higher energy asymmetric "humps" 
(remnants of the strongly suppressed coherence peaks), 
are reproduced here. However, the weaker low-energy 
"shoulders" observed in BSCCO (Refs.0|l are not con- 
sistent with our data. Since the coherence length is much 
smaller in BSCCO than in YBCO, it is more likely that 
mean-field theory will break down in the vortex cores of 
BSCCO, and that strong-correlation effects not captured 
in our model are responsible for the shift of weight to 
higher energies. We do note, however, that the vortex- 
induced states at ±7 meV are seen to persist beyond the 
vortex corCj^ in partial agreement with theory. 



IV. SUMMARY AND CONCLUSIONS 

We find that the negatively charged, antiferromagnetic 
vortex core is preserved up to a long-range Coulomb in- 
teraction strength of V = 0.35 in the underdoped super- 
conductor. For t ^ 0.5 eV, this corresponds to a static 
dielectric constant for the bound charges of et ~ 20. 
Experiments indicate that the static dielectric constant 
of the bound charges is on the order of 10 to 30+ in 
the undoped parent insulatoriS^iS^ and in lightly doped 
cuprates,^® a range which is typical for other oxides. 
This suggests that negatively charged vortices, induced 
by the presence of antiferromagnetic order in the vor- 
tex core, could well be a robust feature of the cuprate 
superconductors in the underdoped region of the phase 
diagram. It is important to note that variation in mate- 



rial properties, such as the chemistry and ordering of the 
interlayer donor regions, or the interlayer distance itself, 
could easily lead to non universal behavior across various 
families of the cuprates. In any case, one of the princi- 
pal results of this paper is that static antiferromagnetic 
order cannot be sustained in the absence of a negatively 
charged core. The nucleation of antiferromagnetic or- 
der as the magnetic field, the doping, or the long-range 
Coulomb interaction strength are changed is always ac- 
companied by a transition to a negatively charged vortex 
core. On this basis, we argue that measurements of the 
vortex core charge could provide strong evidence for — or 
against — antiferromagnetic ordering in the vortex cores 
of high- Tc superconductors. 

Qualitative agreement of the local density of states 
with STM datci2i2i4 is dramatically improved when the 
long-range Coulomb repulsion is included in our model, 
suggesting that the splitting of the zero-bias conduc- 
tance peak [at least in YBCO (Ref. 0)] can be explained 
by the presence of antiferromagnetic order in the vor- 
tex core. We find that the Il-triplet pairing ampli- 
tude is preserved outside of the vortex core, as long 
as there is coexistence of AFM and dSC orders. We 
have also seen that a stripelike ordering of the antifer- 
romagnetism and the charge density can occur when the 
AFM interaction strength is increased, and that this one- 
dimensional anisotropy is strongly enhanced by the long- 
range Coulomb interaction. 

At the relatively high magnetic fields to which we are 
limited by our approach, the vortices are close enough 
to each other that static AFM order, which has a cor- 
relation length ^ ^ ^, is extended throughout the bulk 
of the sample and should be robust against the efl^ects of 
fiuctuations; AFM is still noticeable at the boundary of 
a unit cell size of 26 x 52 sites, which corresponds to a 
field of 19 T. At much lower magnetic fields, where the 
intervortex distance 3> ^, fluctuations may well de- 
stroy the long range AFM order. However, even in the 
absence of long-range order, it is possible that signatures 
of localized antiferromagnetism would be detectable at 
low frequencies in the vicinity of the vortex cores. 

To our knowledge, only one attempt at a direct mea- 
surement of the vortex core charge in high- Tc cuprate su- 
perconductors has been reportedJ£ Nuclear quadrupole 
resonance (NQR) measurements of a field-induced shift 
in the NQR frequency, which is sensitive to the local 
charge density, of the copper nuclei suggest that slightly 
overdoped YBa2Cu307 has negatively charged vortices 
and underdoped YBa2Cu40g has positively charged 
vortices. The results of this experiment suggest a dop- 
ing dependence of the vortex core charge that is opposite 
to that reported in this paper, and in seeming contra- 
diction with the many observations of field-induced anti- 
ferromagnetic orderi^iSiSi^SiSi^SiSi^ in the cuprates and 
with the doping dependence of the Hall signi^ It would 
be interesting to see the results of these and other™ mea- 
surements of the vortex-core charge in the more highly- 
underdoped cuprates. The results of this paper suggest 
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that such measurements could well be definitive. 
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APPENDIX A: CONVERGENCE METHOD 

The long-range Coulomb interaction severely destabi- 
lizes self-consistent calculation of the order parameters, 
by causing divergent long-wavelength oscillations in the 
charge density from one iteration to the next. In order to 
suppress these divergent oscillations we tried a number 
of convergence methods, in increasing order of sophisti- 
cation, which we review below. 

The simplest way to suppress these oscillations is to 



linearly mix the input dxl™'')) the output {\Xout)) 
from the current iteration to generate the input for the 

next iteration |x|r+'^) = (1 " a)lxlr^) + «lxiut^)- While 
this method has been used with some success in other 
self-consistent calculations (see Ref. 'ti', and references 
therein) it is unable to suppress the divergent oscillations 
in the charge density in the BdG calculations of this pa- 
per for any reasonable values of a (the reasons for this 
failure are discussed in Ref. ItS) ■ 

A more sophisticated approach is Broyden's method. 
Let 



(™) \ 
Xout / 



(Al) 



where \F) = |0) at convergence. Expanding \F) to linear 
order near convergence gives 

(A2) 

where the Jacobian 



I. 



(A3) 



If convergence occurs at iteration m + 1, such that 
ji^l^+i)) = 0, then 



(A4) 



where G = [J]-Ms the inverse Jacobian. This equation 
is used to generate a new set of inputs after each iter- 
ation. The goal of Broyden's method is to improve the 
inverse Jacobian at each iteration so that the new inputs 



Ix 



(m+l) 



) are as close to convergence as possible. 



Define, for convenience and for numerical accuracy, 
|^(™)^ _ 



ixir') 
ixir^) 



xir'^)r 

(m-l)^ 



IX 



ixir^) 



IX 



(A5) 
(A6) 



such that (Ax^")|Ax(™)) = 1. At convergence G^™) = 
Q(m-i)^ This suggests that the path to convergence can 
be found by minimizing the change 



E = 



(m-l). 



(A7) 



(where ||...|| denotes the Frobenius norm) in the inverse 
Jacobian from one iteration to the next, subject to the 
constraint [from Eq. HA2|l ] that 



lAx^"') = -G('"^|Ai^(")) 



(A8) 



This minimization leads to a recursion relation for the 
inverse Jacobian (G^™+^^) that is the basis of Broyden's 
method. While the use of Broyden's method does go 
some way towards suppressing oscillations in the charge 
density from iteration to iteration, we found that a long- 
wavelength oscillation (with a period in "time" of two 
iterations) in the charge density would slowly build up 
over many iterations and cause the system to diverge. 

A higher order approach is the modified Broyden's 
method, first introduced by Vanderbilt and Louie (VL)i£^ 
The problem with Broyden's method is that the inverse 
Jacobian G^™' has not been required to satisfy 



(A9) 



for all previous iterations n < m, as it should. 

Following VLS information from all previous itera- 
tions can be introduced by doing a least-squares mini- 
mization of the following "error function:" 



S = L^2||j(™)_j(0)|| 



|Ai^("))-f Jj(")|Ax^"^) 



(AlO) 

where ujn represents the weight given to the information 
from the nth iteration. As an aside, we note that we tried 
various weighting schemes and found that choosing 



UJ„ 



(/p(m)| p(m)\ \ 
\ ' / ] , OJo= 0.01, (All) 
(xlr^ixlrv 



which gives more weight to the information from itera- 
tions in which the difference |xout) ~ Ixin) is small, works 
well for our particular problem. The new error function 
introduced by VL allows the minimization of both the 
change and the error [based on Eq. HA9|) ] in the Jacobian 
from iteration to iteration. Note that we are now consid- 
ering the change in the Jacobian and not the inverse Ja- 
cobian, and that the change in the Jacobian [[J'-™'— Jj'^'H, 
is now defined relative to JJ^"-'. This may seem like a step 
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backwards, since it is the inverse Jacobian that is needed 
to generate the inputs for the next iteration [as in Eq. 
(IA4(l ]. This point was made by Johnson^i^ who intro- 
duced a numericaUy more efficient version of the modified 
Broyden's method by defining an error function in terms 
of the inverse Jacobian G. However, we have found that, 
for our particular problem, using the error function based 
on the Jacobian leads to a procedure which requires far 
fewer iterations to reach self-consistency. Furthermore, 
by making the standard assumptioniiZi that the initial 
Jacobian j'"' = I/a (meaning that linear mixing is used 
to construct the first guess) we are able to use the nu- 
merical innovations of Johnson to greatly enhance the 
efficiency of the method of VL by avoiding any explicit 
calculations of the inverse. 

We start by minimizing the error function [Eq. ljA10|) ] 

with respect to the new Jacobian. Setting dE/dl[^^ = 
gives (as in Ref. ^ 



j(m) ^ p(m)j^(m)j 



where 



and 



3(m) 



in) I 



(A12) 



(A13) 



(A14) 



In order to update the inputs for the next iteration, as 
in Eq. ljA4|l . one needs the inverse Jacobian G*-™^ = 
^(m)[p(m)]-^i^ Assuming that = I/a, 



[r 



("1)1-1 



(A15) 



Borrowing an idea from Johnsonj24 we expand [F^™^]"^ 
to infinite order in the vectors |AF*^"^) and jAx*^")} and 
then resum: 

m 

[r^™^]"' =al + a J2 t^n^On^|AF("))(Ax(')|, (A16) 



where jni = [^I - a]J and a„- = w.Wj (Ax^'^ | AF(j)). 
This procedure exchanges the inversion of a large N x N 
matrix for that of a much smaller m x m matrix. The 
inverse Jacobian is thus 



^^ilki\u'-''^){Ax''^^, (A17) 

k,e=i 



where = Wfe(G(o)|AFW) 

have used the identity 



[Ax'-'^-')) and where we 



UJ, 



, 



to simplify the product /3^™^ [F 

Substituting Eq. HA17|I into the update HA4|I gives 



-6kt H — -Ikt 



(A18) 



fe=i 



where 



,(m) 



Y^lkl 



(A20) 



and 



(A21) 



have been introduced for numerical convenience. The 
modified Broyden's method described here should be gen- 
erally applicable to problems other than the one dis- 
cussed in this paper 

The use of the modified Broyden's method eliminates 
the divergent oscillations in the charge density by tak- 
ing information from all previous iterations to construct 
the input for the next iterationZ&iS. Since we start with 
a converged solution for = and since the largest 
changes at 1^ > are to the charge density we restrict 
application of the modified Broyden's method to updates 
of the Hartree shift. The inputs for all other order pa- 
rameters are taken to be the outputs of the current iter- 
ation. We found that it was necessary to fix the antifer- 
romagnetic order parameter — which strongly infiuences 
the charge density — and allow the Hartree shift to par- 
tially converge over several iterations using the modified 
Broyden's method. Temporarily fixing the antiferromag- 
netism keeps the convergent endpoint of the Hartree shift 
within reach of the modified Broyden's method (which is 
based on the assumption that the starting state is not 
too far from convergence). Once the Hartree shift has 
partially converged, we allow the antiferromagnetism to 
change and then restart the modified Broyden's method 
for the Hartree shift by setting J = I/a and taking as 
initial inputs the most recent outputs. In this way the 
system proceeds stepwise towards a convergent endpoint 
without exceeding the reach of the modified Broyden's 
method in any given step. 
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FIG. 1: From left to right, electron number density relative to the average density 5{ni) = {rn) —0.875, for V = Q, 0.2, and 0.35 
on one half of a 20 x 40 unit cell (the other half is equivalent by symmetry) for an interaction strength J = 1.15. For V = 0.2 
a region of reduced electron density screens the peak in the core. At = 0.35 the vortex core charge and the screening charge 
have changed sign, and modulations in the electron density have a much smaller amplitude (on the order of 10~^ electrons per 
site). The bottom row shows two-dimensional (2D) density plots of the 3D plots in the top row. 
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FIG. 2: From left to right, the AFM order rm (in units of t) for V = 0, 0.2, and 0.35 on one half of a 20 x 40 unit cell, with 
J = 1.15 and riavo = 0.875. With V — 0.2 the magnitude of the AFM has been approximately halved and a,t V — 0.35, where 
the vortex core is no longer negatively charged, the AFM order is negligible. The shape of the AFM order is unaffected by the 
dramatic changes in the structure of the electron density induced by the LRC interaction. The bottom row shows 2D density 
plots of the 3D plots in the top row. 
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FIG. 3; From left to right, the magnitude of the dSC order |Ai| (in units of t) for V = 0, 0.2, and 0.35 on one half of a 20 x 40 
unit cell, with J — 1.15 and riavo = 0.875. The structure of |Ai| is unchanged by the LRC, but the magnitude is gradually 
reduced due to weakening of the nearest-neighbor attraction that generates the dSC pairing. The bottom row shows 2D density 
plots of the 3D plots in the top row. One can see that the coherence length ^ increases as the strength of the d-wave pairing is 
reduced. 




FIG. 4: Peak values of dSC |Ao| (□) far from the vortex 
core, and AFM |mcore| (°) and electron density 5(ncore) (A) 
at the center of the vortex core as a function of increased LRC 
strength V, for a 20 x 40 unit cell with J — 1.15 and riavc ~ 
0.875. All quantities are plotted relative to their V — values. 
Between V = 0.3 and V — 0.35, the density of electrons in 
the core crosses over from positive to negative and the AFM 
order disappears. Lines are a guide to the eye. 
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FIG. 5: Antiferromagnetic order mcorc (dotted lines), in units 
oft, and electron number density (5{ncorc) (solid lines) at the 
center of the vortex core as a function of unit cell size L x 2L 
for LRC strengths of V" = (o), V = 0.15 (□) and V = 0.25 
(0), with J = 1.15 and n^vc = 0.875. When V = the den- 
sity at the core increases toward half filling (0.125 in this plot) 
with increasing unit cell size. For V > the dependence be- 
comes non-monotonic with the peak charge density occurring 
in the 22 x 44 unit cell. The anomalous increase in the AFM 
at V = 0.25 for L — 24, 26 is due to the development of a ID 
anisotropy discussed in Sec. 1111 El Lines are a guide to the 
eye. 




FIG. 6: Antiferromagnetic order mcorc (dotted lines) in units 
of t, and electron number density 5(ncorc) (solid lines) at the 
center of the vortex core as a function of the average den- 
sity riave for LRC strengths of \/ = (o), 1/ = 0.15 (□) and 
V = 0.25 (0) for a 20 X 40 unit cell with J = 1.15. The re- 
gion containing AFM and negatively charged vortices shrinks 
slowly at first and then rapidly toward half filling with in- 
creasing V. Lines are a guide to the eye. 
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FIG. 7: From left to right, electron number density S{ni) — (m) —0.875, AFM order rrn, and dSC order |Ai| for an interaction 
strength J = 1.3 and a LRC strength V = 0.35 on one half of a 20 x 40 unit cell. The bottom row shows 2D density plots of 
the same. Both the electron density and the AFM order show a strong anisotropy, while the dSC order is only mildly affected. 
Note that, in contrast to the J = 1.15 results for V = 0.35, the AFM order survives and the charge density is negative at the 
vortex core. 




FIG. 8: The local density of states at the center of the vortex 
core (top), away from the vortex core along the node direction 
(middle), and averaged over all sites within a coherence length 
from the vortex center (bottom), for J = 1.15 and riavc = 
0.875 on a 26 X 52 unit ceU. Solid lines are for V = 0.15 
and dashed lines are for V = 0. The lines have been shifted 
vertically for clarity. 



